Improving Eeciency of Large Time-scale Molecular Dynamics Simulations of Hydrogen-rich Systems

A systematic analysis is performed on the e ectiveness of removing degrees of freedom from hydrogen atoms and/or increasing hydrogen masses to increase the e ciency of molecular dynamics simulations of hydrogen-rich systems such as proteins in water. In proteins, high-frequency bond-angle vibrations involving hydrogen atoms limit the time step to 3 fs, which is already a factor of 1.5 beyond the commonly used time step of 2 fs. Removing these degrees of freedom from the system by constructing hydrogen atoms as dummy atoms, allows the time step to be increased to 7 fs, a factor of 3.5 compared to 2 fs. Additionally a gain in simulation stability can be achieved by increasing the masses of hydrogen atoms with remaining degrees of freedom from 1 to 4 u. Increasing hydrogen mass without removing the high-frequency degrees of freedom allows the time step to be increased only to 4 fs, a factor of 2 compared to 2 fs. The net gain in e ciency of sampling con gurational space may be up to 15% lower than expected from the increase in time step due to the increase in viscosity and decrease in di usion constant. In principle introducing dummy atoms and increasing hydrogen mass do not in uence thermodynamical properties of the system and dynamical properties are shown to be in uenced only to a moderate degree. Comparing the maximum time step attainable with these methods (7 fs) to the time step of 2 fs that is routinely used in simulations and taking into account the increase in viscosity and decrease in di usion constant, we can say that a net gain in simulation e ciency of a factor of 3 to 3.5 can be achieved. 2 Introduction The maximum time step in Molecular Dynamics (md) simulations is limited by the smallest oscillation period that can be found in the simulated system. Bondstretching vibrations are the fastest atomic motions in a molecule, typically in the order of 10 fs. A classical treatment of these motions is not correct because such vibrations are in their quantum-mechanical ground state. For a proper treatment quantum mechanical calculations should be included. It stands to reason, therefore, that for normal (classical) md these motions are ignored altogether, i.e. they can be better represented by a constraint. For the remaining degrees of freedom, the shortest oscillation period as measured from a simulation is 13 fs for bond-angle vibrations involving hydrogen atoms, see Table I. Taking as a guideline that with a Verlet (leap-frog) integration scheme a minimum of 5 numerical integration steps should be performed per period of a harmonic oscillation in order to integrate it with reasonable accuracy1, the maximum time step will be about 3 fs. This is slightly larger than the 2 fs which is routinely used in md simulations of biomolecules in water. Disregarding these very fast oscillations of period 13 fs (which are also in the quantum-mechanical ground state) the next shortest periods are around 20 fs, which will allow a maximum time step of about 4 fs. In simulations with constrained bond lengths and angles, it has recently been shown that hydrogen atom dihedral angle motions (e.g. rotation of hydroxyl groups) impose a 5 fs limit on the time step, while non-hydrogen atomic collisions (Lennard-Jones \rattling") restrict the time step to a maximum of 10 fs2. The fastest motions in a simulation will invariably involve hydrogen atoms, since these are by far the lightest atoms present in all biological systems. The biologically relevant behavior of these systems mostly takes place on large timescales; at least several nanoseconds but often seconds or even beyond. On these time scales hardly any in uence can be expected from the tens of femtoseconds long oscillations of hydrogen atoms. The obvious solution would be to constrain all bond angles involving hydrogen atoms in all molecules, in addition to all bond lengths. With the shake constraint algorithm3 this can already be done, but shake tends to break down with time steps beyond 2 fs. The generally more robust and faster lincs constraint algorithm4 is now preferred above shake5, but cannot handle the highly connected constraints that arise from constraining both bonds and angles4. The cleanest solution would be to eliminate these high-frequency degrees of freedom from the system altogether. For hydrogen atoms in large molecules (e.g. proteins) this can be implemented in a rather straightforward manner. In stead of connecting a hydrogen atom with bonds, angles and dihedrals to the molecule, the position of the hydrogen will be generated every md step based on the position of three nearby heavy atoms. All forces acting on the hydrogen atom will be redistributed over these heavy atoms. A particle that is treated in this manner is referred to as a dummy atom6, 7. To keep the total mass in the system constant, the mass of each hydrogen atom that is treated in this way, should be added to the bonded heavy atom. Care should be taken that for groups with internal rotation (e.g. hydroxylor amine-groups) only the other internal degrees of freedom of the group should be xed, but the rotational freedom should remain. 3 A special case is the movement of water molecules. The internal geometries of the popular water models are already rigid, so the high-frequency motions are in this case librational motions of the whole water molecule. The only way these motions can be slowed down is by increasing the moments of inertia of the water molecule. This can be done best by increasing the mass of the hydrogen atoms while decreasing the mass of the oxygen atoms, such that the total mass will remain unchanged. A similar modi cation can be made for groups with internal rotational freedom, as hydroxyl or amine groups, which will display motional frequencies close to those of water. All this, of course, constitutes a non-trivial deviation from \physical reality" and requires justi cation. The equilibrium distribution of a system of particles that behaves according to classical statistical mechanics (as in md), is not dependent on the masses of the individual particles. This is not strictly true if constraints are present because of the mass dependence of the metric tensor correction, but these corrections are usually zero or negligible. Therefore we are at liberty to choose appropriate masses without a ecting thermodynamic properties. This idea was originally proposed by Jacucci & Rahman at a c.e.c.a.m. workshop in 19748 and has been investigated in some detail since6, 9, 10 and used in our laboratory before11. None of these studies has been ne grained in the increase of the hydrogen mass. Moreover, those that had a systematical approach, have not kept the total mass of the system constant and have failed to notice that this in fact scales system time. This is most pronounced in Mao et al.12 who report simulations at 600 K and normal time steps with all masses uniformly increased by a factor of 10, but who are actually simulating normal masses with a p10 times smaller time step at 6000 K. This increase in system temperature accounts for all improved conformational sampling that was reported. Also the approach as outlined above has not been applied to simulations of biological (macro)molecules, which require additional modi cations. The issue is, again, the time scales that are of interest for the study of biological systems. On these time scales, at least several nanoseconds, the properties of the water will average and the in uence of water on the dynamics of the system is exerted through its bulk properties, like di usion constant, viscosity and dielectric relaxation time. Hydrogen motions in proteins are almost uncoupled from the main chain vibrations and will therefore hardly in uence the behavior of the system on these time scales. Methods MD Simulations All md simulations were performed using the following parameters and methods, unless stated otherwise. The Verlet integration scheme (leapfrog) was used13. The gromos-87 force eld14 was used, with increased repulsion between water oxygen and carbon atoms15. Explicit hydrogens were de ned for the aromatic rings, the resulting parameter set is the one referred to as SW by Daura et al.16. Periodic boundary conditions with a rectangular box were applied. Lincs4 was used to con4 strain all covalent bonds in non-water molecules. The settle algorithm was used to constrain bond lengths and angles in the water molecules17. The temperature was controlled using weak coupling to a bath18 of 300 K with a time constant of 0.1 ps. Protein and water were independently coupled to the heat bath. Initial velocities were randomly generated from a Maxwell distribution at 300 K, in accordance with the masses that were assigned to the atoms. The pressure was also controlled using weak coupling with a time constant of 1.0 ps. The starting conformation for the water simulations was generated by equilibrating 820 spc (Simple Point Charge) water molecules19 in a 3.2 nm cubic box for 100 ps at 300 K and 1 bar, using a time step of 2 fs and other parameters as stated above. Subsequently, the masses were reassigned (see Table II) and an additional equilibration was done for 10 ps using the same parameters. The same procedure was applied to a smaller (1.9 nm) cubic box containing 216 spc water molecules and an elongated box of 1:9 1:9 5:6 nm containing 648 spc water molecules. The resulting conformations for each hydrogen mass and box shapes and sizes were used as a starting conformation for all simulations of water. Short simulations of 1 ps starting from the equilibrated water structures of the small box of 216 spc waters, for each set of atomic masses were performed to determine the total energy drift as function of atomic mass and time step. Time steps of 0.5 to 15 fs were used; simulations with larger time steps were not stable. A shift function20 was applied for Coulomb and Lennard-Jones interactions which decreases the potential over the whole region and lets potential and force decay smoothly to zero between 0.5 and 0.75 nm. This introduces some artifacts into the simulation20, 21, but it e ectively removes noise from cuto e ects, enabling an accurate assessment of the simulation accuracy as determined by the time step. Neighbor list generation was performed every time step to exclude possible errors from di usion of particles in between neighbor list updates. No pressure coupling was applied, instead the system was equilibrated to the right pressure and density and subsequently simulated at constant volume. No temperature coupling was applied. Longer simulations of 100 ps starting from the equilibrated water structures for each set of atomic masses were performed to determine the dynamical properties of water, i.e. di usion constant, lifetime of hydrogen bonds and viscosity. For determination of the viscosity the elongated box of 648 spc waters was used. For the other determinations the cubic box of 820 spc waters was used. A time step of 2 fs was used. A twin-range cuto for non-bonded interactions was employed with a short-range cuto for Lennard-Jones and Coulomb interactions of 0.9 nm which were calculated every simulation step and a long-range cuto of 1.1 nm for Coulomb interactions which were calculated during neighbor-list generation, every 10 steps (20 fs). Neighbor searching was done based on the centers of geometry of the water molecules22. This is the same parameter set as was used by Van der Spoel et al.23. The simulations of a protein were performed with the small protein HPr (nmr pdb entry 1hdn24). This 85 residue / protein consists of a four-stranded antiparallel -sheet anked on one side by three anti-parallel -helices. The protein was solvated by generating a cubic box of spc water molecules, such that the minimum distance between the protein and the edge of the periodic box would be 0.6 nm, resulting in a cubic box of 4.7 nm. All water molecules from the generated box of water that were within 0.23 nm of a protein atom were removed, leaving 2985 water 5 molecules around the protein. The resulting conformation was energy minimized with harmonic constraints on the atomic coordinates of the protein. Subsequently a round of 10 ps of md was performed, also with harmonic constraints on the atomic coordinates of the protein to relax the water orientation near the protein. The nal conformation was used as starting conformation for simulations of the protein in water. Additionally, another 10 ps of md was performed starting from this conformation. From these two nal conformations the water was removed and the protein was allowed to relax for 1 ps in a vacuum environment. The nal conformations from both 1 ps vacuum simulations were used as starting conformations for simulations of the protein in vacuum. Short simulations of 1 ps of the protein in vacuum were performed to determine the total energy drift as function of time step. Four di erent topology types were used: the normal topology (\normal 1 u"), with the hydrogen atoms fourfold increased in mass (\normal 4 u"), with dummy hydrogens (\dummy 1 u") and with dummy hydrogens and remaining hydrogens fourfold increased in mass (\dummy 4 u"). Time steps of 0.5 to 7 fs were used; simulations with larger time steps were not stable. No cuto for Lennard-Jones or Coulomb interactions was applied and no periodic boundary conditions were used. No temperature coupling was applied. Although simulations of a protein in vacuum are generally not relevant to the majority of applications, these simulations do allow for a relatively accurate estimate of the energy drift, which is not possible for a simulation of the protein in water. Distortions of the shape of the protein by the vacuum environment might in uence energy drift, but are unlikely to occur within the 1 ps duration of the simulations. Long simulations of 1 ns of the protein in water were performed to determine the long-term properties of a protein using the \normal 1 u", \normal 4 u", \dummy 1 u" and \dummy 4 u" topologies, with time steps ranging from 1 to 7 fs. Longrange Coulomb interactions were calculated using pppm25, 26 with a grid spacing of 0.09 nm. Neighbor list generation was performed every 10 time steps. No pressure coupling was applied, instead the system was equilibrated to the right pressure and density and subsequently simulated at constant volume, resulting in average system pressures ranging from 14 to 80 bar with an average of 44 bar. Within this range of pressures, no in uence is likely to exist on the properties of the protein. All md simulations were carried out using the GROMACS molecular dynamics package27, 28 on a Silicon Graphics (sgi) Power Challenge with mips R10000 processors and on sgi O2 Workstations with mips R5000 processors. CPU times for the long (1 ns) runs of the protein in water on the sgi Power Challenge machine are summarized in Table III, for a total of 130 days of cpu time. System Topology Normal topologies, with constraints on all bonds and no constraints on angles, were generated using standard GROMACS topology building tools. These tools were modi ed to optionally produce the modi ed topologies containing the dummy atoms and remove all bond, angle and dihedral de nitions that have become obsolete due to the introduction of the dummy atoms. Also optionally, masses of all remaining normal hydrogen atoms can be increased by a factor of 4, while subtracting this increase from the bonded heavy atom. More details are included in the following 6 section. Construction of Dummy Atoms The goal of de ning hydrogen atoms as dummy atoms is to remove all highfrequency degrees of freedom from them. In some cases not all degrees of freedom of a hydrogen atom should be removed, e.g. in the case of hydroxyl or amine groups the rotational freedom of the hydrogen atom(s) should be preserved. Care should be taken that no unwanted correlations are introduced by the construction of dummy atoms, e.g. bond-angle vibration between the constructing atoms could translate into hydrogen bond-length vibration. Additionally, since dummy atoms are by definition mass-less, in order to preserve total system mass, the mass of each hydrogen atom that is treated as dummy atom should be added to the bonded heavy atom. All forces acting on the dummy atoms must be redistributed over the constructing atoms7. Dummy atom positions can be easily calculated from any three atoms that have a xed orientation with respect to each other. If the constructing atoms move signi cantly with respect to each other, normalized vectors will have to be used to ensure the right position for the dummy atom. This results in complicated derivatives in the force redistribution, which are described in detail in appendix A. Taking into account these considerations, the hydrogen atoms in a protein naturally fall into several categories, each requiring a di erent approach, see also Figure 1: hydroxyl (-OH) or sulfhydryl (-SH) hydrogen: The only internal degree of freedom in a hydroxyl group that can be constrained is the bending of the C-O-H angle. This angle is xed by de ning an additional bond of appropriate length, see Figure 1A. This removes the high frequency angle bending, but leaves the dihedral rotational freedom. The same goes for a sulfhydryl group. Note that in these cases the hydrogen is not treated as a dummy atom. single amine or amide (-NH-) and aromatic hydrogens (-CH-): The position of these hydrogens cannot be constructed from a linear combination of bond vectors, because of the exibility of the angle between the heavy atoms. In stead, the hydrogen atom is positioned at a xed distance from the bonded heavy atom on a line going through the bonded heavy atom and a point on the line through both second bonded atoms, see Figure 1B. planar amine (-NH2) hydrogens: The method used for the single amide hydrogen is not well suited for planar amine groups, because no suitable two heavy atoms can be found to de ne the direction of the hydrogen atoms. In stead, the hydrogen is constructed at a xed distance from the nitrogen atom, with a xed angle to the carbon atom, in the plane de ned by one of the other heavy atoms, see Figure 1C. amine group (umbrella -NH2 or -NH+3 ) hydrogens: Amine hydrogens with rotational freedom cannot be constructed as dummy atoms from the heavy atoms they are connected to, since this would result in loss of the rotational freedom of the amine group. To preserve the rotational freedom while removing the hydrogen bond-angle degrees of freedom, two \dummy masses" are 7 constructed with the same total mass, moment of inertia (for rotation around the C-N bond) and center of mass as the amine group. These dummy masses have no interaction with any other atom, except for the fact that they are connected to the carbon and to each other, resulting in a rigid triangle. From these three particles the positions of the nitrogen and hydrogen atoms are constructed as linear combinations of the two carbon-mass vectors and their outer product, resulting in an amine group with rotational freedom intact, but without other internal degrees of freedom. See Figure 1D. Additionally, all bonds, angles and dihedrals that are de ned on one of the degrees of freedom that were removed, are also removed. This boils down to removing all bonds to dummy atoms, all angles that involve two or three dummy atoms and all dihedrals that involve at least one dummy atom and for which all other atoms are used in constructing the dummy atom(s). Note that this leaves the whole force eld unchanged. As a second option, all remaining hydrogen atoms (i.e. in hydroxyl, sulfhydryl, amine groups and water) can be increased in mass, with the increase subtracted from the bonded heavy atom. This leaves the total mass constant, but increases the moment of inertia of the group, e ectively slowing down the motions. For the dummy mass and dummy atom construction of the amine group (type D as described above) this will have the net result of the dummy masses being placed further apart in accordance with the desired increase in moment of inertia. It should be noted that in the gromos-87 force eld aliphatic hydrogens are implicit, i.e. they are represented as united carbon-hydrogen atoms. For an all-atom force eld with all hydrogen atoms explicit an additional number of hydrogen dummy atoms will have to be constructed every time step, e.g. for the 85 residue protein, 165 hydrogen atoms are present in the gromos-87 force eld but an additional 488 aliphatic hydrogens are implicit. Constructing the 165 hydrogen dummy atoms takes far less than 0.1% of the total computer time, so constructing 653 hydrogen dummy atoms will still have no noticeable e ect on the total cost of simulating. Determination of System Properties Motional Periods Periods of oscillation were measured from a simulation of the protein in water with a time step of 0.5 fs, by taking the highest frequency signi cant peak from the spectrum of the motion. Periods are also calculated from the force eld parameters using: T = 2 I fc 12 (1) where I is the moment of inertia, determined by bond lengths and masses, and fc the force constant of the corresponding angle or dihedral potential. Eq. (1) neglects e ects of coupling with the surroundings. 8 Energy drift As outlined above, the simulations used to determine the drift in the total energy of the simulated system were performed with neither temperature nor pressure coupling. For the water box a shifted potential for electrostatic and Lennard-Jones interactions was used to eliminate cuto e ects and for the protein no cuto for interactions was used. This was done to minimize as much as possible all sources of integration errors (notably cuto e ects) except for those caused by the ratio between time step and fastest motional periods. Also, double precision (8 bytes) oating point calculations were used during the simulation. The limited accuracy of summing up millions of interactions in single precision gives rise to additional drift that obscures the e ects we want to investigate, especially in the smaller time-step regime. As a measure of the accuracy, the root-mean-square (rms) averaged drift in the total energy obtained from a least squares linear t over the last 0.9 ps of 1 ps simulations is used1. Since the drift in the total energy of a well-integrated system is di usive in nature, an appreciable number of independent simulations needs to be performed in order to get an accurate estimate of this drift. The root-mean-square drift provides a reliable measure for the accuracy of the simulation. For the water box, 12 simulations were performed for each time-step/mass combination, each with a di erent random seed to generate initial velocities. For the protein, two di erent starting conformations were used, for which 6 simulations with di erent random seeds were performed for each time-step/topology-type combination. This adds up to 1440 one-ps simulations of the water box, for a total of 12.4 hours cpu time on a mips R10000 processor, and 560 one-ps simulations for the protein, for a total of 5.4 hours of cpu time. The uctuation of the total energy, which can be determined easily from a single simulation, is an inappropriate measure to assess the simulation accuracy of a Verlet-type (leapfrog) integration scheme, since it is of second order in time step, while the drift is of second to third order1. Di usion Constant Di usion constants (D) for water were calculated from the mean square displacement (msd) of the water molecules using the Einstein relation for di usion in three dimensions29: Dkr(t) r(0)k2E = 6Dt (2) D was determined by a linear t to the plot of the msd vs. time. Viscosity The procedure for determination of the viscosity was modi ed after Berendsen30. Viscosity was determined in a non-equilibrium simulation setup where an external shear-stress acceleration eld was applied: ai;x = Acos 2 zi lz (3) 9 with ai;x being the acceleration in the x direction, A the acceleration amplitude, zi the z-coordinate of the particle, lz the length of the box in the z-direction. Application of this shear-stress acceleration gradient induces a velocity gradient of the same shape. For a Newtonian uid the dynamic viscosity is given as the the ratio between the applied acceleration amplitude A and the resulting velocity amplitude V 30: = A V lz 2 2 (4) where is the density and lz the box length in the z-direction of the system. The scaling procedure used in temperature coupling was modi ed to exclude the induced velocity gradient while applying temperature scaling. Care was taken to choose the acceleration amplitude low enough to prevent the appearance of ordering in the water and high enough to get a velocity gradient that is discernible over the thermal velocities. An amplitude of 0.07 nmps 2 was found to perform the best (E. Apol, personal communications); this results in a velocity amplitude of the order of 0.1 nmps 1 which corresponds to roughly 10% of the root-mean-square thermal velocity at 300 K, which is 1.1 nmps 1. For the same reason, the acceleration eld was applied along the longest edge (3 times the length of the other edges) of the rectangular simulation box. The velocity V amplitude was calculated using a spatial Fourier component30: V = 2 N Xi vi;xcos 2 zi lz (5) which was stored at every time step. The viscosity was calculated from eq. (4), using the average velocity amplitude over the last 90 ps of the simulations to exclude the equilibrational part in which buildup of the velocity gradient still occurs. Lifetime of hydrogen bonds Hydrogen bonds between water molecules were de ned using a simple angle and distance criterion, i.e. angle hydrogen-donor-acceptor 60 and distance donoracceptor 0.35 nm, yielding a switch function which is 1 when a hydrogen bond is present and 0 otherwise. The hydrogen-bond lifetime is determined as the half-life time of the autocorrelation of the switch function. Results Water Energy Drift In Figure 2A the energy drift is plotted as a function of hydrogen mass and time step. Note that for clarity the graphs for 2, 3, 6 and 7 fs time step are not shown, they all lie in between the plotted graphs. The drift as a function of time step is of second order, which lies within the expected range of second to third order1. At 10 time steps above 7 fs for the 1 and 8 u simulations and above 10 fs for the 4 and 5 u simulations a transition occurs from second to higher order. Taking a rather arbitrary value of 10 as a maximum for the order of the drift as a function of the time step, a maximum time step can be determined for each hydrogen mass, as is summarized in Table II, these time steps correspond to the sharp increase in the slope of the plots in Figure 2A. For water with the normal mass distribution, the maximum time step is 6.6 fs. The largest attainable time step of 10.4 fs occurs with water with a hydrogen mass of 5 u, an increase of a factor of 1.6. According to eq. (1) this is consistent with the increase in moment of inertia of a factor of 2.5. This correspondence is a clear indication that the librational frequency of (spc) water is the major factor determining the maximum possible time step for accurate integration. Alternatively, a maximum for the magnitude of the drift could be chosen to determine maximum time-steps, however, these will be rather inconsistent due to the large uctuation that is present in the magnitude of the drift. The very sharp increase in the order of the drift allows for a much more accurate determination of maximum time-steps. Di usion The di usion constant as determined from the simulations for di erent atomic masses ranges from 4.1 10 9m2 s 1 at a hydrogen mass of 1 u to 3.3 10 9m2 s 1 at a mass of 4 u, or a decrease of a factor of 1.2 (see Table II). Compared to the di erence between di usion constants of spc water (4.1 10 9m2 s 1) and real water (2.3 10 9m2 s 1), a di erence of a factor of 1.8, the variation caused by changing the hydrogen mass is relatively small. Viscosity Viscosity for di erent atomic masses ranges from 4.3 10 4 kgm 1 s 1 at a hydrogen mass of 1 u to 5.3 10 4 kgm 1 s 1 at a mass of 6 u. At a mass of 4 u, which yields the highest maximum time step, the viscosity is 4.9 10 4 kgm 1 s 1, an increase of a factor of 1.1 (see Table II). Compared to the di erence between the viscosity of spc water (4.3 10 4 kgm 1 s 1) and real water (8.0 10 4 kgm 1 s 1), a factor of 1.9, the variation caused by changing the hydrogen mass is small. Still the result can be signi cant in the sense that large-scale consorted motions in the simulation, e.g. domain motions of proteins, will be limited by the viscosity of water, which means that a higher viscosity of the water will result in slower protein motion. It can be expected that the maximum slowing-down will be of similar order as the increase in the viscosity of the water, i.e. about 14%. Hydrogen bond lifetime The hydrogen bond lifetime increases monotonically with the hydrogen mass, from 0.67 ps at a mass of 1 u to 0.95 ps at 8 u (see Table II), when the single deviating value at 3 u is ignored. At a hydrogen mass of 4 u the lifetime is 0.79 ps, which is an increase of a factor of 1.2 with respect to the value of normal spc water. The lifetime of hydrogen bonds in normal spc water (0.67 ps) may be compared with an 11 experimental estimate of 0.59 ps31, on the basis of uctuations in the anisotropy of molecular polarizability as determined from depolarized Rayleigh scattering measurements. Protein Energy Drift In Table IV a summary is given of the results of the protein simulations, see also Figure 2B. For time steps of 1 fs and below, the drift is di usive in nature, as was the case for the water box, which gives rise to relatively large variations in the determined drift. For larger time steps the drift becomes systematic and is always positive. The order of the drift as a function of time step also lies within the theoretically expected range of second to third order. Surprisingly, the magnitude of the drift is virtually identical for all topology types, the only di erence between them is the time step at which a transition from third to higher order occurs. These transition time-steps are summarized in Table IV, column A. Long-term properties From Table III it is immediately obvious that the average rms deviation from the starting structure in the last 100 ps of each run with respect to the starting structure, the secondary structure content in the last 900 ps and the average total number of hydrogen bonds in the last 900 ps are not noticeably in uenced by the introduction of dummy atoms, heavy hydrogen atoms or large time steps. In Table IV, column B, a summary is given of the maximum time steps for which a 1 ns simulation could be performed without errors. These time steps are somewhat smaller than those obtained from monitoring the energy drift, as summarized in column A. Since the energy drift was determined from simulations in vacuum, it could be expected that this di erence is due to the interactions between protein and water. However, tests of long simulations of the protein in vacuum yield the same maximum time steps as those found for the protein in water. This means that the time-step limiting factors arise from the protein and not from the water, as can also be seen from comparison of the maximum time steps found for the water box (see Table II) and for the protein in vacuum (see Table IV). In Figure 3 inter-protein hydrogen-bond distance and angle distributions are shown. It is clear that the introduction of dummy atoms gives rise to slightly broader distributions. This is to be expected since for normal hydrogens the bond angle can adjust itself to accommodate an optimal hydrogen-bonding conformation. Removing this freedom by constructing the hydrogen as a dummy atom, makes adjustment impossible, giving rise to more sub-optimal hydrogen-bonding conformations and hence a slightly broader distribution. In contrast, no e ect on the distributions is noticeable from changing hydrogen masses or time step. In Figure 4 dihedral angle distributions are plotted for the C-NH+3 dihedral angle, averaged over time steps (Figure 4A) and topology types (Figure 4B). It is clear that neither introduction of dummy atoms in the -NH+3 group, nor increase of 12 mass of the hydrogen atoms nor taking larger time steps has a noticeable e ect on the dihedral angle distributions. Discussion and Conclusion In md simulations of proteins, in which bond lengths are constrained, the usual time step is 2 fs. This is only slightly below the absolute maximum, which is shown to be 3 fs, the 2 fs maximum being a practical limit imposed by the use of the shake algorithm. In order to perform simulations at larger time steps, the hydrogen degrees of freedom should be further restricted. This can be done either by de ning additional constraints, or by treating hydrogen atoms as \dummy" atoms which are constructed from neighboring \real" atoms. We have chosen the latter approach because it (i) avoids problems with constraints in planar groups, (ii) allows the combination of two constrained and one exible angle in a plane (as for the backbone amide proton) and (iii) enables the use of the more stable lincs constraint algorithm in stead of shake to satisfy constraints. The removal of hydrogen degrees of freedom is not expected to cause a noticeable disturbance in the physical behavior of the system on longer time scales because the hydrogen motions are almost uncoupled from the main chain vibrations. This stands in contrast to the strongly coupled heavy-atom bond-angle vibrations that in uence the accessible con gurational space and can therefore not be treated as constraints32. Treating all hydrogens as dummy atoms (\dummy 1 u") allows the time step to be increased by a factor of 2.3 (see Table IV, column B). The bottleneck is now the internal rotation or libration of hydrogen-containing groups and of water molecules (see Table I). The frequencies related to such motions will scale inversely proportional to the square root of the moments of inertia and can thus be decreased by modi cation of the atomic masses. For classical simulations the thermodynamic properties do not depend on the (distribution of) atomic masses. Dynamic properties of a protein on longer time scales will only weakly depend on the mass of hydrogen atoms in the protein and depend on the properties of water through its bulk transport properties. Increase of hydrogen mass by a factor of 4 with simultaneous decrease of the mass of the bonded heavy atom to preserve the total mass of the group (\normal 4 u"), only allows for a modest increase of a factor of 1.3 (see Table IV, column B). Combining the use of dummies with mass modi cation (\dummy 4 u") allows for an increase in time step of 2.3, which is identical to that observed for \dummy 1 u". It appears that no additional gain comes from increasing hydrogen masses in a system where most hydrogen atoms in the protein are already treated as dummy atoms, however, a gain in simulation stability is to be expected. The viscosity of water increases and di usion coe cient decreases by roughly 15%; therefore the net gain in simulation e ciency can be up to 15% lower than expected based on the increased time step. An additional factor can also result in an e ciency gain that will be slightly less; when a neighbor list is used, which is to be updated for example every 20 fs, counting in integration steps it must be updated more frequently. Using both dummy atoms and modi ed masses, the next bottleneck is likely 13 to be formed by the improper dihedrals (which are used to preserve planarity or chirality of molecular groups) and the peptide dihedrals. Although the improper dihedrals could be replaced by dummy atom constructions or their potential function be modi ed to reduce the resonance frequencies, the peptide dihedral cannot be changed without a ecting the physical behavior of the protein. Thus we have approached the limit of what can be achieved without a ecting the physical behavior. We would like to conclude from this discussion that measuring the drift in total energy of a simulation allows one to determine the maximum time-step given a maximum order of the drift as a function of the time step. Table IV, column A, shows that this criterion would allow for a time step of 3 fs for normal simulation, and an increase of a factor of 2 for simulating with hydrogen atoms of 4 u. Introducing dummy hydrogens atoms will allow for a gain in maximum time-step of factor of 2.7, irrespective of the mass of remaining hydrogen atoms. As is evident when comparing columns A and B in Table IV, in real-life examples the maximum timesteps will be somewhat lower. It appears that monitoring the total energy drift fails to capture some important features of the simulated system which determine the integration accuracy and stability. This is most pronounced for the \normal 4 u" case; based on energy drift a maximum time step of 6 fs would be expected, but an actual simulation of a protein in water remains stable only up to time steps of 4 fs. It seems therefore best to choose an important property (or a number of properties) of the system for which a reference value or distribution is known (either from experiment or from an accurately performed reference simulation), and monitor this property during simulation. This, however, gives rise to a new problem because for many systems, such a property might be hard to nd. The most practical approach to determine the maximum time-step is simply to determine for which time steps the simulation itself will remain stable. From Table IV it can be seen that increasing hydrogen atom mass will allow for a modest increase in time step from 3 to 4 fs, introducing dummy hydrogen atoms, however, allows the time step to be increased to 7 fs and the combination of increased hydrogen atom mass and dummy hydrogen atoms will give the same time step of 7 fs, but with a slightly less uctuations in various simulation parameters, and presumably a better long-term stability. In a concluding summary we can say that an increase in time step from 2 to 7 fs, a factor of 3.5, for routine md simulations of proteins in water can be achieved by constructing hydrogen atoms in the protein as dummy atoms, leading to a gain in simulation e ciency of a factor of 3 to 3.5. Additional simulation stability can be gained by increasing the mass of all hydrogen atoms with remaining degrees of freedom from 1 to 4 u. Acknowledgments We would like to thank Dr. M. L. C. E. Kouwijzer for stimulating discussions and to A. K. Mazur for some very helpful suggestions. Thanks also go to our undergraduate students for conducting a number of preliminary studies: M. Bergsma, R. M. de Blouw, S. J. Boot, A. Duursma, M. van Faassen, S. Harkema, D. R. Hekstra, M. Mol 14 and C. Smit. B. H. acknowledges support from the Netherlands Foundation for Chemical Research (SON) and K. A. F. acknowledges support from the Foundation for Life Sciences (SLW), both with nancial aid from the Netherlands Organization for Scienti c Research (NWO). 15 A Redistribution of Forces on Dummy AtomsDummy atoms are virtual particles with position rd, which are constructed fromthe positions of the real particles ri. Therefore every rd is a known function ofri's. Any force F d on a dummy atom is redistributed to the real atoms on whichrd depends. When a linear combination of three atoms is used in constructing thedummy atom, the weight for redistributing the forces are equal to those used inthe linear combination. This redistribution becomes nontrivial if normalization isused in constructing the dummy atom, as is the case for the dummy types used foraromatic, amide and amine hydrogens (see Figure 1B and C).The force acting on atom i (F0i) as a result of the force on the dummy atommust be calculated from the partial derivative of the position of the dummy atomwith respect to the position of atom i7:F0ix = @rd@xi @V@rd = @rd@xi F d(6)Here V is the potential energy expressed in positions of real and dummy atompositions. Analogous expressions are valid for the y and z component.For type B (see Figure 1B), the position of the dummy atom d is calculatedfrom the position of the constructing atoms i, j and k as follows:rd = ri + b rij + arjkjrij + arjk j(7)where rij = rj ri. Using eq. (6) to calculate the redistributed force for atoms i,j and k yields:F0i = F d (F d F 1)F0j = (1 a) (F d F 1)F0k =a (F d F 1) where=bjrij + arjk jF 1 = rid F drid rid rid(8)For type C (see Figure 1C) the position is calculated using:rd = ri + b cos rijjrij j + b sin r?jr?j where r? = rjk rij rjkrij rij rij (9)with corresponding forces, using eq. (6):F0i = F db cosjrij j F 1 + b sinjr?j rij rjkrij rij F 2 + F 3F0j =b cosjrij j F 1b sinjr?j F 2 + rij rjkrij rij F 2 + F 3F0k =b sinjr?j F 2where F 1 = F d rij F drij rij rij , F 2 = F 1 r? F dr? r? r? , F 3 = rij F drij rij r?and r? as de ned in eq. (9)(10)16 References1 A. K. Mazur, J. Comp. Phys. 136, 354 (1997).2 A. K. Mazur, J. Phys. Chem. 102, 473 (1998).3 J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comp. Phys. 23, 327(1977).4 B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comp.Chem. 18, 1463 (1997).5 H. J. C. Berendsen, Molecular Dynamics Simulations: The Limits and Beyond,in Computational Molecular Dynamics: Challenges, Methods, Ideas, edited byP. Deu hard, J. Hermans, B. Leimkuhler, A. Mark, S. Reich, and R. D. Skeel,Springer-Verlag, 1998, ISBN 3-540-63242-5.6 C. M. Bennet, J. Comp. Phys. 19, 267 (1975).7 H. J. C. Berendsen and W. F. van Gunsteren, Molecular Dynamics Simulations:Techniques and Approaches, in Molecular Liquids: Dynamics and Interactions,edited by A. J. Barnes, W. J. Orville-Thomas, and J. Yarwood, NATO ASIC 135, pp. 475{500, Reidel, Dordrecht, The Netherlands, 1984, ISBN: 90-277-1817-2.8 G. Jacucci and A. Rahman, The Possibility of Using a Larger Timestep in Molec-ular Dynamics Studies of Water, in Report on Workshop Methods in MolecularDynamics Long Timescale Events, pp. 32{40, Orsay, 1974, C.E.C.A.M.9 D. W. Wood, Computer Simulation of Water and Aqueous Solutions, in Water:A Comprehensive Treatise, edited by F. Franks, volume 6, p. 279, Plenum Press,New York, 1979.10 R. Pom es and J. McCammon, Chem. Phys. Lett. 166, 425 (1990).11 E. Egberts, S. Jan Marrink, and H. J. C. Berendsen, Eur. Biophys. J. 22, 423(1994).12 B. Mao, G. Maggiora, and K. C. Chou, Biopolymers 31, 1077 (1991).13 L. Verlet., Phys. Rev. 34, 1311 (1967).14 W. F. van Gunsteren and H. J. C. Berendsen, Gromos-87 manual, Biomos BV,Nijenborgh 4, 9747 AG Groningen, The Netherlands, 1987.15 A. R. van Buuren, S. J. Marrink, and H. J. C. Berendsen, J. Phys. Chem. 97,9206 (1993).16 X. Daura, B. Oliva, E. Querol, F. X. Avil es, and O. Tapia, PROTEINS: Struct.Funct. Gen. 25, 89 (1996).17 S. Miyamoto and P. A. Kollman, J. Comp. Chem. 13, 952 (1992).18 H. J. C. Berendsen, J. P. M. Postma, A. DiNola, and J. R. Haak, J. Chem.Phys. 81, 3684 (1984).17 19 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans,InteractionModels for Water in Relation to Protein Hydration, in IntermolecularForces, edited by B. Pullman, pp. 331{342, D. Reidel Publishing Company,Dordrecht, 1981.20 K. F. Lau, H. E. Alpher, T. S. Thacher, and T. R. Stouch, J. Phys. Chem. 98,8785 (1994).21 P. E. Smith and B. M. Pettitt, J. Chem. Phys. 95, 8430 (1991).22 R. R. Gabdoulline and C. Zheng, J. Chem. Phys. 16, 1428 (1996).23 D. van der Spoel, P. J. van Maaren, and H. J. Berendsen, J. Chem. Phys. 108,10220 (1998).24 N. A. J. van Nuland, I. W. Hangyi, R. C. van Schaik, H. J. C. Berendsen, W. F.van Gunsteren, R. M. Scheek, and G. T. Robillard, J. Mol. Biol. 237, 544(1994).25 R. W. Hockney and J. W. Eastwood, Computer simulation using particles,McGraw-Hill, New York, 1981.26 B. A. Luty, M. E. Davies, I. G. Tironi, and W. F. van Gunsteren, Mol. Sim.14, 11 (1994).27 H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comp. Phys. Comm.91, 43 (1995).28 D. van der Spoel, A. R. van Buuren, E. Apol, P. J. Meulenho , D. P. Tiele-man, A. L. T. M. Sijbers, R. van Drunen, and H. J. C. Berendsen, gromacsUser Manual version 1.5, Nijenborgh 4, 9747 AG Groningen, The Netherlands.Internet: http://rugmd0.chem.rug.nl/~gmx, 1997.29 M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, OxfordScience Publications, Oxford, 1987.30 H. J. C. Berendsen, Transport properties computed by linear response throughweak coupling to a bath, in Computer Simulations in Material Science, editedby M. Meyer and V. Pontikis, pp. 139{155, Kluwer, 1991.31 C. J. Montrose, J. Bucaro, J. Marshall-Coakley, and T. Litovitz, J. Chem. Phys.60, 5025 (1974).32 W. F. van Gunsteren and H. J. C. Berendsen, Mol. Phys. 34, 1311 (1977).33 D. R. Lide, editor, CRC Handbook of chemistry and physics : a ready-referencebook of chemical and physical data, CRC Press, Boca Raton [etc.], 72 edition,1991.34 W. Kabsch and C. Sander, Biopolymers 22, 2577 (1983).18 List of TablesI Characteristic oscillation periods of atomic motions in md simula-tions. fc: force constant; I : moment of inertia, or atomic mass forbond stretching; calc.: calculated from eq. (1); sim.: highest fre-quency signi cant peak in spectrum of angle respectively dihedralmotion from simulation. An entry of \ " means not applicable, ornot determinable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21II Atomic masses in water; corresponding smallest moments of inertiaI ; resulting dynamical properties: viscosity ; di usion constant D,values of H2O and D2O from Lide et al.33 and hydrogen-bond lifetimeH bond, value of H2O from Montrose31; rms drift of the total energyover 12 runs at a time step of 4 fs; maximum time step ( tmax) ata maximum order of the drift as a function of time step of 10. . . . . 21III Summary of long (1 ns) simulations of the protein in water for simu-lations with \normal 1 u", \normal 4 u", \dummy 1 u" and \dummy4 u" topologies. Simulation parameters: time step; number of steps;total runtime on an sgi Power Challenge with mips R10000 pro-cessors, averaged over the four topology types. Long term averageproperties: rms deviation of all backbone atoms with respect to thestarting structure, averaged over last 100 ps; secondary structurecontent (% of residues not in random-coil conformation, according tothe dssp program34) averaged over 100 to 1000 ps; number of inter-protein hydrogen bonds averaged over 100 to 1000 ps. Entries of \ "indicate failure of the simulation to run without errors, this was alsothe case for time steps larger than 7 fs. . . . . . . . . . . . . . . . . . 22IV Summary of maximum time steps tmax for which A: the drift of thetotal energy as a function of the time step in short (1 ps) simulationsis still of third order; B: long (1 ns) simulations can be performedwithout errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2219 List of Figures1 Schematic view of the di erent types of dummy atom constructionsused. The atoms used in the construction of the dummy atom(s) aredepicted as black circles, dummy atoms as grey ones. Hydrogens aresmaller than heavy atoms. A: xed bond angle, note that here thehydrogen is not a dummy atom; B: in the plane of three atoms, withxed distance; C: in the plane of three atoms, with xed angle anddistance; D: construction for amine groups (-NH2 or-NH+3 ), see textfor details. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 rms averaged drift in the total energy as a function of time step ofA: a box of spc water (648 atoms) with di erent hydrogen masses(graphs for 2, 3, 6 and 7 u are omitted for clarity). B: a small protein(805 atoms) with di erent topology types (\normal 1 u", \normal4 u", \dummy 1 u" and \dummy 4 u"). . . . . . . . . . . . . . . . . 233 Inter-protein hydrogen-bonds: A donor-acceptor distance distribu-tion and B hydrogen-donor-acceptor angle distribution, averaged overall simulations without dummy atoms (\normal") and with dummyatoms (\dummy"). All distributions appear to be insensitive tochanges in time step and hydrogen masses. . . . . . . . . . . . . . . 234 Dihedral angle distributions inLys24-NH+3 of the protein. A: aver-aged over all time steps for each topology type; B: averaged overtopology types for time steps of 1, 2, 4 and 7 fs (see Table III for anoverview of simulations at di erent time steps and topology type). . 2420 Table I: Characteristic oscillation periods of atomic motions in md simulations. fc: forceconstant; I: moment of inertia, or atomic mass for bond stretching; calc.: calculated fromeq. (1); sim.: highest frequency signi cant peak in spectrum of angle respectively dihedralmotion from simulation. An entry of \ " means not applicable, or not determinable.fcIPeriod (fs)motion:(kJmol 1) (u nm2) calc. sim.bond stretch, H400 000m = 1 u 1010bond stretch, heavy atoms 500 000m = 12 u 3020water libration0:005928water rotation0:00591300angle, H3750:0103220angle, heavy atoms4500:2715445angle-NH+3 group, C-N-H3750:0103222angle-NH+3 group, H-N-H7500:0102313improper, planar16728improper, tetrahedrical33527dihedral, peptide-bond330:2048928dihedral,-NH+3 group3:80:023 48989dihedral, OH group1:30:0094 5343Table II: Atomic masses in water; corresponding smallest moments of inertia I; resultingdynamical properties: viscosity ; di usion constant D, values of H2O and D2O from Lideet al.33 and hydrogen-bond lifetime H bond, value of H2O from Montrose31; rms drift ofthe total energy over 12 runs at a time step of 4 fs; maximum time step ( tmax) at amaximum order of the drift as a function of time step of 10.mass (u)IDH bond drift Etottmax (fs)H O (u nm2) (10 4 kgm 1 s 1) (10 9m2 s 1)(ps)kJmol 1 ps 1(fs)1 160:00594:34:080:671:046:62 140:01044:73:890:740:868:93 120:01334:93:790:890:4210:04 100:01484:93:340:790:3610:3580:01485:13:500:840:4710:4660:01335:33:350:840:598:6740:01045:23:340:880:437:5820:00595:13:600:950:615:6real H2O8:02:30:59real D2O2:021 Table III: Summary of long (1 ns) simulations of the protein in water for simulationswith \normal 1 u", \normal 4 u", \dummy 1 u" and \dummy 4 u" topologies. Simulationparameters: time step; number of steps; total runtime on an sgi Power Challenge with mipsR10 000 processors, averaged over the four topology types. Long term average properties:rms deviation of all backbone atoms with respect to the starting structure, averaged overlast 100 ps; secondary structure content (% of residues not in random-coil conformation,according to the dssp program34) averaged over 100 to 1000 ps; number of inter-proteinhydrogen bonds averaged over 100 to 1000 ps. Entries of \ " indicate failure of thesimulation to run without errors, this was also the case for time steps larger than 7 fs.t Nsteps cpurms deviation (nm)sec. struct. (%)# H-bonds(fs) ( 103) timenormaldummynormal dummynormaldummy(hours) 1 u 4 u 1 u 4 u 1 u 4 u 1 u 4 u 1 u 4 u 1 u 4 u1 1000341 0:11 0:12 0:18 0:11 87 86 85 84 114 112 117 1142 500171 0:15 0:15 0:17 0:13 86 85 87 85 114 111 116 1183 333114 0:15 0:09 0:17 0:13 87 87 82 86 114 115 114 1164 250860:16 0:15 0:1487 85 86114 116 1145 200680:12 0:1687 85116 1146 167570:13 0:1588 85117 1187 143490:11 0:1889 84113 116Table IV: Summary of maximum time steps tmax for which A: the drift of the totalenergy as a function of the time step in short (1 ps) simulations is still of third order; B:long (1 ns) simulations can be performed without errors.topology type tmax (fs)A Bnormal 1 u33normal 4 u64dummy 1 u87dummy 4 u87 DBbbα ACFigure 1: Schematic view of the di erent types of dummy atom constructions used. Theatoms used in the construction of the dummy atom(s) are depicted as black circles, dummyatoms as grey ones. Hydrogens are smaller than heavy atoms. A: xed bond angle, notethat here the hydrogen is not a dummy atom; B: in the plane of three atoms, with xeddistance; C: in the plane of three atoms, with xed angle and distance; D: construction foramine groups (-NH2 or-NH+3 ), see text for details. 22 11010−210−1100101102103 DriftEtot(kJmol−1ps−1 )Water1 u4 u5 u8 uy=x2y=x3 110Proteinnormal 1 unormal 4 udummy 1 udummy 4 uy=x3 Timestep (ps)Figure 2: rms averaged drift in the total energy as a function of time step of A: a boxof spc water (648 atoms) with di erent hydrogen masses (graphs for 2, 3, 6 and 7 u areomitted for clarity). B: a small protein (805 atoms) with di erent topology types (\normal1 u", \normal 4 u", \dummy 1 u" and \dummy 4 u"). 0.250.30.35Distance (nm)0204060Angle (O)normaldummy ABFigure 3: Inter-protein hydrogen-bonds: A donor-acceptor distance distribution and Bhydrogen-donor-acceptor angle distribution, averaged over all simulations without dummyatoms (\normal") and with dummy atoms (\dummy"). All distributions appear to beinsensitive to changes in time step and hydrogen masses.23 −120 0120Topology Typenormal 1 unormal 4 udummy 1 udummy 4 u −120 0120Timestep1 fs2 fs4 fs7 fs Dihedral Angle (O)ABFigure 4: Dihedral angle distributions inLys24-NH+3 of the protein. A: averaged over alltime steps for each topology type; B: averaged over topology types for time steps of 1, 2, 4and 7 fs (see Table III for an overview of simulations at di erent time steps and topologytype).24